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ABSTRACT 

High energy materials such as polymer bonded explosives are commonly used as propellants. These 
particulate composites contain explosive crystals suspended in a rubbery binder. However, the explosive 
nature of these materials limits the determination of their mechanical properties by experimental means. 
Therefore micromechanics-based methods for the determination of the effective thermoelastic properties 
of polymer bonded explosives are investigated in this research. Polymer bonded explosives are two- 
component particulate composites with high volume fractions of particles (volume fraction > 90%) and 
high modulus contrast (ratio of Young's modulus of particles to binder of 5,000-10,000). Experimentally 
determined elastic moduli of one such material, PBX 9501, are used to validate the micromechanics 
methods examined in this research. The literature on micromechanics is reviewed; rigorous bounds on 
effective elastic properties and analytical methods for determining effective properties are investigated in 
the context of PBX 9501. Since detailed numerical simulations of PBXs are computationally expensive, 
simple numerical homogenization techniques have been sought. Two such techniques explored in this 
research are the Generalized Method of Cells and the Recursive Cell Method. Effective properties 
calculated using these methods have been compared with finite element analyses and experimental data. 
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INTRODUCTION 

Large scale simulations of the mechanical response of containers filled with high energy 
(HE) materials require knowledge of the mechanical properties of the materials involved. 
Though mechanical properties of HE materials can be determined experimentally, hazards as- 
sociated with experiments on these materials, as well as the attending costs, make this option 
unattractive. As computational capabilities have grown and improved numerical techniques 
developed, numerical determination of the properties of HE materials has become possible. In 
this research, we explore micromechanics-based numerical methods for the determination of 
the mechanical properties of HE materials, specifically polymer bonded explosives (PBXs). 
The polymer bonded explosive PBX 9501 has been characterized experimentally and thus pro- 
vides a basis for validating numerical calculations. 

PBXs provide unique challenges for micromechanical modeling. These materials are vis- 
coelastic particulate composites, contain high volume fractions of particles, and the modulus 
contrast between the particles and the binder is extremely high. For example, PBX 9501 con- 
tains about 92% by volume of particles and the modulus contrast between particles and binder 
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of 20,000 at room temperature and low strain rates. In this research it is assumed that PBXs are 
two-component particulate composites with the particles perfectly bonded to the binder. The 
components of PBXs are assumed to be isotropic and linear elastic, and only the prediction 
of elastic moduli and coefficients of thermal expansion (CTEs) of PBXs is addressed. Finite 
element (FEM) analyses are used to compute the effective properties of PBX 9501 models and 
compared with properties computed using the Generalized Method of Cells (GMC). A new 
technique called the Recursive Cell Method (RCM) has been developed to address the limi- 
tations of GMC. Effective properties of PBX 9501 are also calculated using RCM for some 
microstructures simulating PBX 9501. 

POLYMER BONDED EXPLOSIVES 

Polymer bonded explosives are mixtures containing high volume fractions of stiff explo- 



sive crystals suspended in a continuous compliant binder. Some common PBXs (Gibbs and 
Popolato 19801 are shown in Table [T] 



TABLE 1. Common polymer bonded explosives and their components. 



PBX 


Crystal 


Weight 


Volume 


Binder 


Weight 


Volume 






(%) 


(%) 




(%) 


(%) 


PBX 9501 


HMX 


95 


92 


Estane 5703 + BDNPA/F 


5 


8 


PBX 9010 


RDX 


90 


87 


KEL-F-3700 


10 


13 


PBX 9502 


TATB 


95 


90 


KEL-F-800 


5 


10 



PBX 9501 

PBX 9501 is a mixture of HMX particles coated by a binder composed of Estane and a 
plasticizer (BDNPA/F). The HMX crystals are in the stable (3 phase, have a monoclinic struc- 
ture (Skidmore et al. 1997), and are linearly elastic at room temperature. The binder in PBX 
9501 is viscoelastic. However at room temperature and low strain rates, the response of the 
binder is close to linear elastic. The composite (PBX 9501) is also linear elastic at room tem- 
perature and low strain rates. The elastic properties of PBX 9501 and it's components at 23°C 
and a strain rate of 0.05/s ( |Gray III et al. 1998j|Zaug 1998[ ) are shown in Table |2j 



TABLE 2. Elastic properties of PBX 9501 , HMX and binder. 



Material 


Young's Modulus 


Poisson's Ratio 


CTE 




(GPa) 




(xl0~ 5 /K) 


PBX 9501 


1.03 


0.35 


12 


HMX 


15.3 


0.32 


11.6 


Binder 


0.001 


0.49 


20 



MICROMECHANICS METHODS 

The term "micromechanics" describes a class of methods that use continuum mechanics 
for determining the effective material properties of composites given the material properties 
of the constituents. Methods of determining the effective elastic properties and coefficients of 
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thermal expansion of particulate composites, containing up to 70% by volume of particles, are 
well established (Has hin 1983[ Markov 2000} Buryachenko 2001 1. Rigorous bounds on effec- 
tive properties ( [Milton 2002 1 and various analytical methods for calculating effective proper- 
ties ( |Torquato 2001 ) are also available. However, these techniques are not particularly useful 
for high volume fraction (> 90%) and high modulus contrast composites such as PBX 9501. 
Therefore numerical techniques are essential to the solution of micromechanics problems for 
polymer bonded explosives. 

Rigorous Bounds and Analytical Methods 

Hashin-Shtrikman (Hashin and Shtrikman 1963 ) and third order (Milton 198 1 ) bounds have 



been calculated for the effective elastic moduli of PBX 9501 . The Rosen-Hashin bounds ( Rosen 



and Hashin 1970 1 on the effective coefficient of thermal expansion of PBX 9501 have also been 



calculated. These bounds are shown in Table [3] The bounds on the elastic moduli are quite far 
apart and hence of no practical use. However, the bounds on the CTE are within 1% of each 
other. 



TABLE 3. Bounds on the effective properties of PBX 9501 . 





Expt. 


Hashin-Shtrikman Bounds 


Third Order Bounds 


Upper 


Lower 


Upper 


Lower 


Bulk Modulus (GPa) 


1.1 


11.4 


0.15 


11.3 


0.22 


Shear Modulus (GPa) 


0.4 


5.3 


0.01 


5.0 


0.07 


CTE(xlO" 5 /K) 


12 


12.3 


11.6 







Analytical methods of interest for high volume fraction particulate composites include the 



composite spheres assemblage (Hashin 1962 1, the self-consistent scheme (SCS) (Berryman 



and Berge 1996), and the differential effective medium (DEM) approach ( Markov 2000[ ). Each 



of these methods makes simplifying assumptions about the microstructure of the composites. 
The CTE can be calculated given knowledge of the isotropic bulk modulus of the compos- 
ite (Rosen and Hashin 1970). The effective elastic moduli of PBX 9501 calculated using SCS 
and DEM, and the corresponding CTEs are shown in Table [4j It can be observed that the ana- 
lytical methods do not predict effective elastic moduli that are close to the experimental values. 
However, the predicted effective CTE is close enough to the experimental value to be of use 
and no further numerical calculation is required for this property. 



TABLE 4. Effective properties of PBX 9501 from analytical models. 





Bulk Modulus 


Shear Modulus 


CTE 




(GPa) 


(GPa) 


(xlO" 5 /K) 


Experiment 


1.1 


0.4 


12 


Self-Consistent Scheme 


11.0 


4.7 


12.9 


Differential Effective Medium 


0.2 


0.08 


12.5 



Numerical Approximations 

The effective elastic moduli of a composite can be determined approximately by solving 
the governing differential equations using numerical methods. This process involves the de- 
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termination of a Representative Volume Element (RVE), the choice of appropriate boundary 
conditions and the solution of the resulting boundary value problem. The effective stiffness 
tensor {C*- kl ) of the composite is then calculated from the relation 

r au dV = C* ikl [ e kl dV (1) 



V JV 

where V is the volume of the RVE, <7jj are the stresses, and e^i are the strains. 

Finite element analysis is the most commonly used numerical technique use to determine 
effective properties. The difficulties involved in discretization and solution of actual particulate 
geometries has led to simulations of simple geometries such as square or hexagonal particles, 
mostly in two dimensions. With improvement in computational capabilities random distribu- 
tions of particles are being simulated more frequently (Mishnaevsky and Schmauder 2001 



The large computational cost associated with finite element analyses and the poor accuracy in 
areas of high stress gradient have led to the exploration of alternative methods of calculating 



effective properties. Discrete spring network models (Day et al. 1992 1, integral equation based 



methods (Greengard and Helsing 1998) and Fourier transform based methods (Michel et al 
|2001) show the most promise. However, a large amount of computation is still required for 
convergence in all these methods. 

Generalized Method of Cells ( GMC) 

The Generalized Method of Cells (GMC) ( |Aboudi 1996] ) has been used to model the mi- 



cromechanical behavior of various types of composites with relative success. The advantage of 
this method over other numerical techniques is that the full set of effective elastic properties can 
be calculated in one step. In this technique, the RVE is discretized into a number of subcells as 
shown in Figure [T] Continuity of displacements between subcells and subcell equilibrium are 
satisfied in an average sense using integrals over subcell boundaries. GMC has been shown to 
be more computationally efficient than finite elements for modeling fiber composites. A refor- 
mulated version of GMC (Pindera and Bednarcyk 1999 1 has been used for the calculations in 
this research. 



Y,2 

p Subcell 




FIG. 1 . Schematic of the Generalized Method of Cells. 
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Recursive Cell Method (RCM) 

The GMC technique underestimates the effective properties of PBX 9501. The Recursive 
Cell Method (RCM) has been developed to provide a computationally efficient and accurate 
alternative to GMC. A schematic of RCM is shown in Figure [2] The RVE is discretized into 
subcells as in GMC. However, instead of calculating effective properties of the whole RVE 
in a single step, the effective properties of small blocks of subcells are determined at a time. 
The effective properties of the RVE are calculated by combining the effective properties of 
blocks using a recursive process. The effective properties of each block of subcells may be 
determined using any accurate numerical technique. We use a finite element based technique 
in this research. The RCM recursive scheme has been found to reduce the computational cost 
and remedy the shear-coupling problem of GMC. 



RVE - Level 



RVE - Level 1 




a 
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Subcell 



HiBBEIEr- 

I I I 
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Block 




Final RVE RVE - Level 2 

FIG. 2. Schematic of the Recursive Cell Method. 



SIMULATION OF PBX 9501 MICROSTRUCTURE 

Random sequential packing techniques were used to generate five different microstructures 
based on the particle size distribution in the dry blend of PBX 9501 (Skidmore et al. 19971. 
Each simulated microstructure contains around 100 particles that occupy approximately 87% 
of the volume of the square RVE, as shown in Figure [3] Since particles occupy 92% of the 
volume in PBX 9501, the binder is assumed to be "dirty", i.e., it contains 40% HMX particles 
by volume so that the total volume fraction of HMX in the composite is 92%. The effective 
properties of the dirty binder were calculated using the differential effective medium (DEM) 
approach. 

The two-dimensional plane strain effective stress-strain relation for the composite is given 

by 

(2) 
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FIG. 3. Sample simulated PBX 9501 microstructure. 



The effective stiffness matrix terms (C?) of the five PBX 9501 microstructures have been 
calculated using finite elements (FEM), GMC, and RCM. The average effective properties 
for the five microstructures and the corresponding standard deviations are shown in Table [5] 
On average, the five microstructures exhibit square symmetry. The FEM calculations (with 
256 x 256 square elements) overestimate the effective normal stiffness (C^ and C| 2 ) of PBX 
9501 by around 15% to 25% on average. The corresponding effective stiffness matrix terms 
obtained by RCM (with 8x8 blocks, each containing 32 x 32 subcells) are around 60% higher 
than the FEM values and around two times the experimental values for PBX 9501. On the 
other hand, the values obtained from GMC form the normal stiffness terms are less than l/10th 
the FEM and the experimental values. The FEM calculations underestimate the values of C\ 2 
for PBX 9501 by around 40%. RCM calculates the same average value of this component 
of the stiffness matrix as FEM. However, GMC predicts a much lower value that is around 
l/7th of the stiffness of PBX 9501. The shear stiffness component (C* m ) of PBX 9501 is 
calculated quite accurately by FEM and the error is only around 10%. RCM underestimates 
this stiffness component by around 60%. However, the value predicted by GMC is around 
1/300 the experimental value and quite close to the Hashin-Shtrikman lower bound on the 
shear modulus. 



TABLE 5. Effective stiffness of simulated PBX 9501 microstructures. 







Std. Dev. 




Std. Dev. 


Cl2 


Std. Dev. 


°66 


Std. Dev. 




(GPa) 


(GPa) 


(GPa) 


(GPa) 


(GPa) 


(GPa) 


(GPa) 


(GPa) 


PBX 9501 


1.60 


1.60 


0.88 


0.38 


FEM 


2.03 


0.81 


1.85 


0.51 


0.54 


0.20 


0.51 


0.34 


RCM 


3.20 


0.96 


2.99 


0.50 


0.54 


0.27 


0.27 


0.16 


GMC 


0.15 


0.02 


0.15 


0.02 


0.12 


0.01 


0.01 


0.001 



DISCUSSION 

The FEM calculations on the simulated PBX 9501 microstructures show that the mi- 
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crostructures are nearly isotropic. A better approximation of isotropy can be obtained if larger 
RVEs are modeled. The variability between models can also be reduced with larger RVEs. 
However, larger RVEs require larger meshes to model the increased degrees of freedom and 
hence are limited by computational cost. The FEM calculations are also found to predict values 
of stiffness that are slightly higher than the experimentally determined stiffness of PBX 9501. 
Three-dimensional calculations could lead to the prediction of lower stiffness for the PBX 9501 
microstructures. In addition, the actual material contains voids and cracks that reduce the stiff- 
ness. The finite element calculations could incorporate interface elements to model interfacial 
debonds. 

The RCM calculations can be performed at least three times faster than the full finite el- 
ement calculations. Higher computational speeds can be obtained by reducing the number of 
subcells per block at each level of recursion though there is some loss in accuracy. Though 
the normal stiffness matrix components predicted by RCM are higher than those predicted 
by FEM, the values are quite acceptable considering the large modulus contrast between the 
components of PBX 9501. Additionally, the RCM calculations provide much better estimates 
of the effective properties than GMC, analytical models or rigorous bounds at relatively low 
computational cost. 

The low effective normal stiffness predicted by GMC is because stress bridging effects 
are underestimated by this technique. The low shear stiffness is due to the lack of coupling 
between the normal and shear stresses and strains in GMC. Modifications of GMC that attempt 
to rectify these problems have been found to lose the computational efficiency of the original 
formulation. 

CONCLUSION 

The effective properties of the models of PBX 9501 are quite close to experimentally de- 
termined properties of PBX 9501. Therefore model RVEs using circular particles distributed 
according the particle distribution of a dry blend of PBX 9501 can be use to simulate PBX 
9501. The RCM technique can be used to calculate effective properties of particulate com- 
posites quite accurately and at low computational cost. The technique is much faster than full 
FEM calculations using the same amount of discretization. Improved performance and accu- 
racy can be obtained in RCM by proper choice of the number of subcells per block at each 
level of recursion. GMC is inadequate for calculating the effective properties of high energy 
materials because of the underestimation of stress bridging effects and the lack of coupling 
between normal and shear stresses and strains. 
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